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ABSTRACT 

A flat distribution of low gas entropy in the core region of galaxy clusters is a feature com- 
monly found in Eulerian cosmological simulations, at variance with most standard simulations 
of Smoothed Particles Hydrodynamics fashion. From the literature it is still unclear whether 
this difference is entirely due to numerical artifacts (e.g. spurious transfer from gravitational 
energy to thermal energy), to physical mechanisms (e.g. enhanced mixing in Eulerian codes) 
or to a mixture of both. This issue is related to many still open lines of research in the charac- 
terization of the dynamical evolution of the baryons in galaxy clusters: the origin of the cool 
core/non-cool core bi-modality, the diffusion of metals within galaxy clusters, the interplay 
between Active Galactic Nuclei and the Intra Cluster Medium, etc. 

In this work, we aim at constraining at which extent the entropy core is affected by nu- 
merical effects, and which are the physical reasons for its production in cosmological runs. 
To this end, we run a set of 30 high resolution re- simulations of a ~ 3 • 10 14 M o //i clus- 
ter of galaxies with a quiet dynamical history, using modified versions of the cosmological 
Adaptive-Mesh-Refinement code ENZO and investigating many possible (physical and nu- 
merical) details involved in the production of entropy in simulated galaxy clusters. 

We report that the occurrence of a flat entropy core in the innermost region of mas- 
sive cluster is mainly due to hydro-dynamical processes resolved by the numerical code (e.g. 
shocks and mixing motions) and that additional spurious effects of numerical origin (e.g. ar- 
tificial heating due to softening effects) affect the size and level of the entropy core only in a 
minor way. 

Using Lagrangian tracers we show that the entropy profile of non-radiative simulations 
is produced by a mechanism of "sorting in entropy " which takes place with regularity during 
the cluster evolution. The evolution of tracers illustrates that the flat entropy core is caused 
by physical mixing of subsonic motions (mostly driven by accreted sub-clumps) within the 
shallow inner cluster potential. 

Several re- simulations were also produced for the same cluster object with the addition 
of radiative cooling, uniform pre-heating at high redshift (z = 10) and late (z < 1) thermal 
energy feedback from AGN activity in the cluster, in order to assess the effects of such mecha- 
nisms on the final entropy profile of the cluster. We report on the infeasibility of balancing the 
catastrophic cooling (and recovering a flat entropy profile) by means of the investigated trials 
for AGN activity alone, while for a sub-set of pre-heating models, or AGN feedback plus pre- 
heating models, a flat entropy distribution similar to non-radiative runs can be obtained with 
a viable energy requirement. Complementary analysis are presented also for a major merger 
cluster, obtaining similar results and achieving a generally good consistency with X-ray data 
for the entropy distribuion in real galaxy clusters. 
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1 INTRODUCTION 



plete absence in the core of galaxy clusters simulated with La- 
grangian approaches (such as Smoothed-Particle-Hydrodynamics 
codes, SPH), has been object of an interesting debate in the last 
few years and among a number of different groups (e.g. Frenk et 
al.1999; Voit et al.2005; Dolag et al.2005; Lin et al.2006; Wadsley 
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et al.2008; Tasker et al.2008; Mitchell et al.2008; Springel 2010; 
Abel 2010). 

The issue of the real inner entropy profile of non-radiative 
galaxy clusters may be loosely connected to the case of real clus- 
ters, which are interested by cooling and energy feedback by sev- 
eral astrophysical sources; nonetheless, this topic represents one of 
the main diagnostic to compare cosmological simulations of galaxy 
clusters performed with different numerical approaches. 

The evolution of gas in the simulated Universe, when only 
shock heating is present as explicit source of energy dissipation, 
represents a rather simple scenario to study the thermodynamics 
of cosmic baryons starting from simple initial conditions. Under- 
standing the systematics that affect the generation of entropy in the 
different codes, even in this rather idealized setup, would represent 
a useful step forward in all research topics dealing with the numer- 
ical treatment of non-reversible processes in astrophysical plasmas 
at all scales. In non radiative simulations of galaxy clusters, the 
task of understanding the correct entropy distribution in the inner- 
most cluster regions is made complicated by a few circumstances: 
a) the radius of the entropy core produced in grid codes is rather 
small r CO re ~ 0.1R V i r (where R V i r is the cluster virial radius), 
and even the most resolved cosmological simulations can concen- 
trate only a moderate number of resolution elements, of the order 
of N ~ 10 2 - 10 4 , inside this region: therefore resolution and sam- 
pling problems may always be present; b) the baryon accreted by 
a cluster are interested by several dynamical processes across their 
evolution (shock heating, violent relaxation, gas mixing and slosh- 
ing of the Dark Matter peak) which distribute entropy in clusters 
in different ways; disentangling the various effects within the same 
galaxy cluster is usually not a trivial task; c) due to the typical ra- 
dial entropy distribution in galaxy clusters, physical mixing driven 
by matter sub-clumps in-falling from the outside regions, and nu- 
merical mixing or spurious heating from N-body noise would have 
the same net effect, leading to an increase of gas entropy within 
the densest regions in clusters; unluckily, the different numerical 
methods are prone to numerical mixing in a ways difficult to quan- 
tify. Therefore very similar entropy configurations may be degen- 
erate respect to various interplays between physical and numerical 
effects along the whole cluster evolution, and specific numerical 
tests aiming at the close comparison between re- simulations of the 
same objects with different numerical methods are highly desirable 
in this respect. 

From the literature, the first clear indication of a fundamental 
difference between the results of SPH codes and grid ones in galaxy 
cluster simulations was presented in the Santa Barbara Compari- 
son project (Frenk et al.1999). In this work, evidence was obtained 
that the innermost entropy radial distribution in cluster simulated 
in grid methods such as ENZO contain a nearly isoentropic core 
inside ~ 0.1R V i r , at variance with SPH codes; later works basi- 
cally confirmed this trend also at higher resolutions (e.g. Voit et 
al.2005; Wadsley et al.2008). Several reasons were suggested to in- 
terpret this discrepancy: over-mixing in grid codes (e.g. Wadsley et 
al.2008), spurious N-body heating from DM particles in the clus- 
ter core (e.g. Lin et al.2006; Springel 2010), Galilean invariance 
in the gravity solver of grid codes (Tasker et al.2008; Robertson 
et al.2010), lack of physical mixing in SPH codes (e.g. Dolag et 
al.2005; Agertz et al.2007; Abel 2010), pre-shocking in SPH (e.g. 
O'Shea et al.2005) etc. To date, the most detailed analysis of the 
generation of entropy has been presented by Mitchell et al.(2009), 
by studying idealized cluster mergers with a non-cosmological 
setup with the SPH code GADGET2 (Springel et al.2005) and the 



grid code FLASH (Fryxell et al.2000). The authors provided strik- 
ing evidences that the mechanism at work in setting the different 
entropy level between SPH and grid codes takes place at the time 
of the closest encounter between the colliding structures, and it is 
related to the suppression of mixing in SPH because of artificial 
viscosity, which highly suppresses hydro instabilities and mixing 
motions respect to grid codes. It would be now interesting to extend 
the results of this seminal paper to fully cosmological simulations, 
and to the case of clusters in which the mass is assembled with 
realistic and different dynamical history (major mergers or regular 
smooth accretions). 

The present paper is first devoted to constrain at which ex- 
tent all the numerical effects cited above may affect also the spe- 
cific distribution of gas entropy in realistic and high resolution cos- 
mological simulations. Secondly, this paper is devoted to focus on 
the physical mechanisms for the generation and the spreading of 
gas entropy within clusters, in relation to gravitational mechanisms 
(e.g. shock heating and mixing motions) and to non-gravitational 
ones (e.g. radiative cooling and energy feedbacks from astrophys- 
ical sources such as AGNs). To this end we produced a set of 
30 re- simulations of the same galaxy cluster, with a final mass of 
M w 3.1 • 1O 14 M //i and a very quiet dynamical history for most 
of its evolution (complementary results for a major merger cluster 
are shown in the Appendix). 

The first part of this work (Sec(3j explores many of the pos- 
sible numerical mechanisms which may lead to the formation of 
an inner entropy core in cosmological cluster simulations, by using 
customized re- simulations with the adaptive mesh refinement code 
ENZO (Norman et al.2007). In detail, SecEUdiscusses the role of 
the mesh refinement strategy adopted in the simulation; Sec l3.3l es- 
timates the role played by cold unresolved flows and N-body grav- 
itational heating, making changes to the "dual energy formalism" 
method adopted in ENZO; Sec l3.4| investigates the influence of the 
softening length in the computation of the gravitational force in 
the PM method; Sec l3. 5 1 compares re- simulations adopting differ- 
ent maximum resolution to compute hydro-dynamics of baryons. 

In Sec l4. li the physical generation of entropy and the volume 
spreading of it during cluster evolution is investigated by means of 
Lagrangian tracers, injected and evolved in non radiative runs. 

Radiative cooling can completely alter the above picture, in 
systems characterized by a cooling time shorter compared to the 
cluster lifetime (e.g. Katz & White 1993; Fabian 1994). Indeed the 
hot gas phase in the core region of these systems is removed by the 
radiative losses, and the inward motion of the cooling gas would 
produce the theoretical "cooling flow" scenario. However, drastic 
cooling flows are not observed in real clusters, and additional non- 
gravitational heating mechanisms are need to restore (or keep) the 
cooling gas on an higher adiabat (e.g. Evrard & Henry 1991; Kaiser 
1991; Lloyd-Davies, Ponman & Cannon 2000). Two of the most 
promising scenarios in this respect are the pre-heating scenario 
(e.g. White 1991), in which gas is heated before it collapses within 
structures, or AGN feedback (e.g. Churazov et al.2000), in which 
the inner entropy is raised by the energy released, through different 
channels, out of the regions surrounding the growing super massive 
black hole. Many works investigated the above mechanisms with 
cosmological simulations, by implementing non-gravitational heat- 
ing mechanisms in SPH or in grid methods (e.g. Borgani et al.2002; 
Borgani et al.2005; Heinz et al.2006; Sijacki & Springel 2006; 
Younger & Bryan 2007; Burns et al.2008; Mc Carthy et al.2009; 
Teyssier et al.2010). However comparing the outcomes of these 
studies is made complex by the great number of assumption and 
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parameters often involved in the modelization of feedback sources. 
Also considering the underlying fundamental differences of the 
hydro methods discussed above, the characterization of the inner 
entropy profile of radiative cluster with non-gravitational heating 
mechanisms at play is a still very open topic of research for simu- 
lations. For instance, qualitatively similar models applied to GAD- 
GET (Borgani et al.2002) or ENZO simulations (Younger & Bryan 
2007) has lead to different conclusion about the efficiency of pre- 
heating at the scale of galaxy clusters or galaxy groups. The abil- 
ity of cluster merger in quenching (or slowing down) the cooling 
catastrophe in radiative simulations is also a debated issue (Burns 
et al.2008; Poole et al.2008). 

In the second part of the paper (SecHJ) we focus on the addi- 
tional physical mechanisms which are able to affect the gas entropy 
distribution in real galaxy cluster, by studying in detail the effects 
of a) radiative cooling (Sec l4.2l ); b) non-gravitational heating by a 
"uniform" heating mechanism in the early Universe (Sec l4.2.TT ); c) 
non-gravitational heating by a central AGN with outflows, within 
an already formed cluster (Sec j4.2.2l4.2.31 >. Our goal is not that of 
constraining the most likely extra heating mechanism at work in 
real galaxy cluster, but rather to show what is the net effect of plau- 
sible heating mechanism on the entropy distribution of a cluster 
with an ongoing cooling flow, using a budget for the energy release 
under control. 

Section [5] finally summarizes our discussions of the results 
and our conclusions, while in the Appendix complementary tests 
studying the numerical and physical generation of entropy in major 
merger cluster are reported for completeness. 



2 NUMERICAL CODE AND SETUP 

The computations presented in this work were performed us- 
ing the ENZO code, developed by the Laboratory for Computa- 
tional Astrophysics at the University of California in San Diego 
(http://lca.ucsd.edu). 

ENZO is an adaptive mesh refinement (AMR) cosmologi- 
cal hybrid code highly optimized for high performance computing 
(Norman et al.2007 and references therein). It uses a Particle-Mesh 
method to follow the dynamics of the collision-less Dark Matter 
(DM) component (Hockney & Eastwood 1981), and a Eulerian 
solver based on the Piecewise Parabolic Method (PPM, Woodward 
& Colella, 1984). 

The adopted simulational setup is the same as in Vazza, 
Gheller & Brunetti (2010). In summary, cosmological initial con- 
ditions were produced with nested grid/DM particle distributions 
of increasing in order to achieve high DM mass resolution in the 
region of cluster formation; also an implemented mesh refinement 
scheme was applied to trigger mesh refinement based on gas/DM 
over-density and/or velocity jumps across cells. In this scheme, a 
normalized 1-D velocity jump across 3 adiacent cells in the scan 
direction (at a given refinement level) is recursively computed as 
as S = \Av/vmin\, where |v m i n | is the minimum velocity, in ab- 
solute value, among the 3 cells. The scheme is made manifestely 
non Galileian invariant by the presence of Vm%n\ this problem is 
unavoidable in this kind of simulations, because in principle every 
forming shock wave moves on a different reference frame, and a 
run-time procedure to account for this would represent a too large 
computational effort. In order to have this effect under control we 
performed many convergence tests with idealized and cosmological 



Table 1. Main characteristics of the performed runs. Column 1: run identi- 
fication code; column 2: cell resolution at the maximum refinement level; 
column 3: maximum softening length; column 4: mesh refinement strat- 
egy:"D"= gas/DM over-density refinement; "V"=velocity jumps refine- 
ment; column 5: additional numerical parameter of the simulations. r\ \ , 772 
and M t hr are parameters involved in the dual energy formalism switch 
(Sec l3.3h "PH" means that uniform pre-heating has been adopted, assum- 
ing an entropy increase of So at z = 10; "J" means that e je t extra energy 
has been injected at z ^ 1 within the cluster, assuming AGN feedback. 



ID Max Res. [kpc/h] e soft . [kpc/h] AMR note 



da 
KU 




jU 


U V 


non-radiative 


Rl 


25 


25 


DV 


non-radiative 


R2 


25 


50 


D 


non-radiative 


R3 


25 


25 


D 


non-radiative 


R4 


12.5 


25 


DV 


non-radiative 


R5 


12.5 


25 


D 


non-radiative 


R6 


12.5 


12.5 


DV 


non-radiative 


R7 


12.5 


12.5 


DV 


non-radiative 


R18 


25 


100 


D 


non-radiative 


R19 


25 


100 


DV 


non-radiative 


R20 


25 


12.5 


D 


non-radiative 


R21 


50 


50 


D 


non-radiative 


R22 


50 


50 


DV 


non-radiative 


R8 


25 


25 


DV 


771 = 10- 2 


R9 


25 


25 


DV 


771 = 10- 4 


R10 


25 


25 


DV 


m = 1 


Rll 


25 


25 


DV 


772 = 10- 2 


R12 


25 


25 


DV 


M thr = 1.1 


R13 


25 


25 


D 


M thr = 1.5 


R15 


25 


25 


DV 


cooling 


R16 


25 


25 


D 


cooling 


PHI 


25 


25 


DV 


cool.+PH(10/ceVcra 2 ) 


PH2 


25 


25 


DV 


cool.+PH(100/ceVcm 2 ) 


PH3 


25 


25 


D 


cool.+PHQOO/ceVcm 2 ) 


PH4 


25 


25 


DV 


cool.+PH(200fceVcm 2 ) 


Bl 


25 


25 


DV 


cool.+J(10 58 ergs) 


B2 


25 


25 


DV 


cool.+J(10 59 ergs) 


B3 


25 


25 


D 


cool.+JQO 59 ergs) 


B4 


12.5 


12.5 


DV 


cool.+J(10 58 er#s) 


B5 


25 


25 


DV 


cool.+PH+J(2 • 10 57 ergs) 



simulations with ENZO, finding showing that a very good numeri- 
cal convergence is achieved in cluster simulations by fixing 8 = 3 
(more detailed discussions can be found in Vazza et al.2009; Vazza, 
Gheller & Brunetti 2010). 

All the re- simulations presented here focus on the evolution 
of the same galaxy cluster, which was the most massive one of the 
sample already presented in Vazza, Gheller & Brunetti (2010). 

In order to role of various numerical effects (such as res- 
olution, softening in the gravitational force, etc.) various re- 
simulations starting with the same initial conditions were produced. 
Table |2 lists the details of all runs performed for this project. The 
re- simulations employing additional physics (such as pre-heating 
and AGN feedback) where produced with original implementations 
made starting from the public 1.5 version of ENZO. 
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Figure 3. Radial profiles for the run R0, showing gas density and DM density profiles (left panel) and gas temperature and gas entropy (right panel). In the 
right panel, the solid lines refer to mass-weighted profiles, while the dashed lines refer to volume-weighted profiles. The vertical arrows show the softening 
length adopted in this run. 
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Figure 1. Redshift evolution of the Dark matter plus gas mass (black), and 
of total gas mass (blue) inside the virial radius of the galaxy cluster studied 
in this work. 



3 NUMERICAL EFFECTS ON THE ENTROPY 
PRODUCTION. 

3.1 The fiducial run 

The fiducial cluster run of this project is R0, which has the same 
numerical setup and mesh refinement strategy already adopted in 
other works by our group (Vazza et al.2009; Vazza, Gheller & 
Brunetti 2010; Vazza et al. 2010). The assumed cosmology is 
the "Concordance" ACDM model, with parameters: Qo = 1.0, 
Qbm = 0.0441, Qdm = 0.2139, Ov = 0.742, Hubble parame- 
ter h = 0.72 and a normalization for the primordial density power 
spectrum ag = 0.8. The mesh refinement on gas/DM over-density 
is triggered since the beginning of the simulation inside a cubic 



with the side of « 5R vir « 12Mpc/h\H z = 30, while the ad- 
ditional mesh refinement triggered by velocity jumps is activated 
from z — 2. 

The lowest resolution level (I = 0) inside the AMR region 
is 220kpc/h, while the maximum refinement level (/ = 3) is set 
to 2bkpc/h. The DM matter particles have a mass resolution of 
6.7 • 1O 8 M //i, and the comoving gravitational softening is e — 
50kpc/h. 

This fiducial cluster run does not include radiative cooling or 
any heating mechanism other than gravitational collapse or shock 
heating. 

Figfl] shows the redshift evolution of the total mass and of the 
gas mass inside the virial cluster region, measured with a spheri- 
cal over-density method. The trend of the cluster mass growth and 
the visual inspection of movies of its evolution confirms that this 
cluster does not experience any violent merger event for z < 1.5, 
and roughly ~ 70 per cent of its mass has been already assem- 
bled at z ~ 0.6. At z = the cluster has a total mass of 
« 3.1 • 1O 14 M /h , an average temperature of 2.2keV and a virial 
radius of Rv = lMMpc/h. 

In FigOwe show the maps of gas density and gas entropy 
through the center of the cluster. The cluster is quite regular in 
shape, with sharp circular shock structures and a well defined inner 
entropy floor of size r cor e ~ lOOkpc/h surrounded by the steep in- 
crease of the ICM entropy, up to a maximum at about ~ 2Mpc/h 
from the cluster center. This cluster may be considered as a "pro- 
totype" of relaxed clusters in the local Universe with intermediate 
mass, as produced by non-radiative cosmological simulations of 
Eulerian fashion. For complementary results on larger mass clus- 



1 In the following, we will refer to this region as to the 'AMR region'. 

2 All throughout this paper, we will refer to S = as to the "gas en- 
tropy", as usually done in cosmological numerical simulations, where P is 
the gas pressure, p is the gas density within a cell and 7 is the adiabatic 
index. 
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Figure 2. Maps of gas density (top) and gas entropy (bottom) for a slice 
crossing the center of the cluster at z = (run RO). The side of the image 
is « h.bMpc/h and the width of the slice is one cell = 25kpc/h. 



Figure 4. Slices showing the absolute value of the 3-D velocity field for 
run RO (Top panel) and run R2 (Bottom panel) at z = 0. The size of the 
image is as in Fig [2] 



ters simulated with the same approach, we refer the reader to recent 
high-resolution re- simulations presented in Vazza et al.(2010), and 
references therein. 

Figure [5] shows the radial profiles of gas density, DM den- 
sity, temperature and entropy for run RO at z = 0. In the right 
panel of this Figure we compare gas mass weighted profiles and 
volume weighted profiles for the cluster run: the two estimates pro- 



vide consistent results within < 2Mpc/h from the cluster center. 
Unless specified, in what follows we will make use of gas density 
weighted averages. 

We investigated the "robustness" of this cluster representation 
by adopting several changes in this simple setup, which will be 
discussed in the next Sections. 
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Figure 5. Radial profiles of kinetic energy for run RO (red solid) and for run 
R2 (blue solid), and of "turbulent" kinetic energy for the same runs (dashed 
lines). The 4 top lines refer to the total kinetic/turbulent energy inside a 
given radius, while the 4 bottom lines refer to the values within shells of 
width 25kpc/h at the same radii. The vertical arrow show the softening 
length adopted in these runs. 



3.2 The role of the mesh refinement strategy. 

Shock heating during the gravitational collapse in the forming clus- 
ter has a leading role in the production of the baseline entropy dis- 
tribution within clusters, since no other physical mechanisms of ir- 
reversible heating are present in simple non-radiative cosmological 
simulations (e.g. Voit et al.2005). Previous works in the literature 
(Dolag et al.2005; Iapichino & Niemeyer 2008; Vazza et al.2009; 
Vazza, Gheller & Brunetti 2010) have suggested that an increase of 
the level of chaotic motions within simulated clusters is expected 
as soon as the effects of a coarse resolution or artificial viscosity 
are limited by ad-hoc techniques, and as a result typically higher 
entropy level is found in the core region of cosmic structures. Also, 
the different ability in modeling shocks and mixing in the central 
phase of cluster mergers has proved to be responsible for the bulk 
of the difference between SPH or grid simulations of galaxy clus- 
ters (Mitchell et al.2009). Obtaining a good spatial resolution of 
chaotic and mixing motion in the ICM is therefore crucial and the 
mesh refinement strategy outlined in Sec l3.ll is designed for this 
purpose (see also Iapichino & Niemeyer 2008). 

We compare the cluster entropy profile in the fiducial run (R0) 
against a re- simulation using the standard mesh refinement based 
on gas/DM over-density alone (run R2). Figure [4] shows the maps 
of the absolute value of the velocity field in the two runs, for a slice 
crossing the center of the cluster at z = 0. Similarly to what re- 
ported in Vazza et al.(2009), the extra refinement strategy allows 
us to preserve a very accurate description of the low density outer 
accretion region, and on the mixing motions following the crossing 
of satellites even if the involved over-density is not enough to trig- 
ger other refinements. The radial distributions of the kinetic energy 
and "turbulent" energy (we fiducially consider as "turbulent" the 
components of the 3-D velocity field which are characterized by 
a coherence scale smaller than < 200kpc/h, see Vazza et al.2009 
Sec.4.2 for further details) are shown in Fig[5] As expected, R0 
shows a slightly higher level of kinetic energy all across the cluster 



Figure 6. Gas entropy profile for run R0 (red) and run R2 (blue). The ad- 
ditional black lines show the fits profiles presented in Voit et al.(2005), for 
S oc r 11 (dashed) and S oc r 12 (dot-dashed). The vertical arrow show 
the softening length adopted in these runs. 

volume, and an kinetic energy in small-scale chaotic motions by a 
factor ^2 — 3 compared to the standard run R2. 

The entropy profiles of the two simulations are shown in Fig(6] 
The additional black lines show the fit profile presented in Voit 
et al.(2005), where a systematic study of entropy profile in non- 
radiative cluster simulations (of SPH and grid fashion) were pre- 
sented. The run with the additional refinement strategy produces a 
significantly larger entropy within the core region, by a ~ 20 — 30 
per cent (the gap is 30 — bOkeVcm 2 ). It is interesting to notice that 
the gas entropy and the gas kinetic/turbulent energy of run R0 is 
larger compared to run R2 within approximately the same radius, 
r < 200kpc/h (the softening for the gravitational force here is 
50kpc/h at the maximum refinement level). 

Already from this, one would speculate that enhanced turbu- 
lent mixing motions are responsible for the level of gas entropy in- 
side clusters; we will focus on this issue with more detail in Sec|4] 



3.3 The role of gravitational N-body heating and cold 
unresolved flows. 

It has been suggested in the recent past that the entropy core in 
grid simulation may be due to an incorrect handling of the gas 
thermal energy/entropy, in the case of poorly resolved cold flows 
dominated by kinetic energy (e.g. supersonic flows in the rarefied 
Universe) and due to noise-induced heating in virialised structures 
(Lin et al.2006; McCarthy et al.2007; Springel 2010). This second 
mechanism may arise because the small-scale velocity fluctuations 
induced by the N-body gravitational field can be readily dissipated 
by the mesh-based hydrodynamics, causing a spurious heating of 
the gas. In Springel (2010) the initial conditions for the Santa Bar- 
bara Comparison (Frenk et al.1999) were evolved with the AREPO 
code using an energy-entropy switch scheme in order to decide if 
the gas entropy had to be updated as the difference between the total 
energy and the kinetic energy (as customarily done in grid codes) 
or by directly using the entropy equation instead. This was done in 
order to ensure that in the case of transonic or subsonic motions 
(M < M thresh , with Mthresh = 1.1) the evolution of gas entropy 
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Figure 7. Comparison of the entropy profiles for different choices of the 771 (R8,R9, left panel) and for different choices of 772 (R10,R11, right panel). 
Additional green profiles are shown for ENZO re- simulations which adopted an additional switch condition on the Mach number to trigger the dual energy 
formalism (R12,R13). The vertical arrow in both panel shows the softening length adopted in these runs. 



was not dominated by any spurious dissipation, since a conserva- 
tive entropy equation was followed. The level of core entropy in 
the simulated cluster was found to be considerably reduced after 
the use of this switch, and the conclusion was that the entropy in 
the core region of clusters simulated with standard methods in grid 
simulations may be significantly affected by this numerical effect. 

This issue is very interesting and deserves investigation; to this 
end we re- simulated run R0 with different choices for the com- 
putation of gas thermal energy. The strategy customairly adopted 
in ENZO is to adopt a "dual energy formalism" (Ryu et al.1993; 
Bryan et al.1995) designed to compute the evolution of gas thermal 
energy whenever the internal energy, e, cannot be calculated cor- 
rectly as the difference between total, E, and kinetic energy (due to 
the fact that E > > e in highly supersonic flows, and round-off er- 
rors may dominate). We briefly recall here its design: a first switch 
condition is followed to update the gas pressure in the case of very 
cold flows 

• p = p(>y-l).(E- v 2 /2) if (E - v 2 )/E > 771 

• p = p(>y-l)-eif(E-v 2 )/E<ri 1 

where p is the gas pressure, v is the modulus of the velocity 
field and 7 = 5/3. The second switch is adopted to update the 
internal gas energy without advecting numerical errors from each 
cell's local neighborhood: 

• e=(E- v 2 /2) if p(E - v 2 )/max[pE] neigh > 772 

• e = v/p{l - 1) tfp(E - v 2 )/max[pE] neigh < 772, 

where max\pE\ neig h is the maximum total energy in the (1- 
D) neighborhood of the cell. This approach ensures that e is not 
contaminated by errors advected by the total energy formulation 
(Bryan et al.1995). However we note that in the above switch con- 
ditions the presence of v makes them non Galilelian invariant, be- 
cause bulk velocities affect the exact value of E and of max [pE] . 
To keep this problem under control, one need to resort to conver- 
gence studies by varying the numerical parameters involved in the 



switch; the customary values set in ENZO cosmological simula- 
tions are 771 = 10 -3 and 772 = 0.1. 

As a first step to investigate the role played by the dual en- 
ergy formalism in the production of the entropy level in clusters, 
we performed several tests by re- simulating the fiducial run adopt- 
ing different choices of the threshold values involved in the above 
"switches": 771 = 10" 2 (R8), 771 = 10" 4 (R9), or 772 = 1 (R10), 
772 = 10~ 2 (Rll). 

The motivation for this kind of tests is that if the low entropy 
gas sitting the in the cluster center at z = is reminiscent of the en- 
tropy production prior to the cluster virialization, then any artifact 
present in the "cold" (T < 10 5 K) Universe at early redshifts (i.e. 
due to an incorrect handling of the gas internal energy, or to features 
related to non Galileian invariance) will be highlighted by different 
choices of the entropy switch. It should be noted that, in absence 
of a re-heating UV background due to stars/AGN (e.g. Haardt & 
Madau 1996), the baryon temperature outside structures can be as 
low as T rsj 1 — 10 K in simulation of this kind. 

By comparing the results of these runs with the fiducial one, 
we report that the net effect of the above changes is quite small, 
as shown in Fig[7] the maximum difference is found for run Rll 
(772 = 10 -2 ) and results only in a ^ 10 per cent larger entropy 
inside the cluster core. Cold unresolved flows are unlikely to be the 
reason of the entropy core in grid simulations. 

As a second step to investigate further the possible role played 
by weak shocks and extra-heating related to N-body gravitational 
noise, we supplemented the dual energy formalism in ENZO with 
an additional Galilelian-invariant switch based on the Mach num- 
ber, with a procedure similar to Springel (2010). 

As shown above the standard switch condition of the dual en- 
ergy formalism may be affected by the presence of large bulk ve- 
locities, which enter in the total energy E, and thus the value of 
the parameter 772 cannot be readily realted to the real Mach number 
of the flow. For this reason we implemented an on-the-fly shock in 
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ENZO simulations, which allows us to apply a Galileian invariant 
switch condition (based on the Mach number) to decide whether or 
not the total energy equation must be used to update the gas dynam- 
ics. This method relies on the temperature jumps method of Ryu et 
al.(2003) and adopts the following procedure: a) candidate shocked 
cells are identified according to the V • v < criterion; b) local gra- 
dients of gas entropy, \7S and gas temperature VT are evaluated, 
and a shock condition is matched whenever VT • VS > 0; c) the 
shock Mach number is estimated by inverting the maximum tem- 
perature jump across the shocked cell: 



softening 



T 2 /Ti 



(5M 2 - 1) • (M 2 + 3) 
16 • M 2 " 



(1) 



This procedure is followed on the fly during the simulation, 
and in this method the total energy equation is used to compute 
the gas internal energy only when p(E — v 2 )/max[pE] ne igh < 
r]2 and M < M thresh simultaneously. This procedure forces the 
code to update the gas thermal energy only in an adiabatic way, 
for flows with M < M t hr', however we note that this formulation 
is not totally equivalent to the energy-entropy formalism discussed 
in Springel (2010), because the internal energy equation does not 
explicitly conserve gas entropy (however it conserves gas energy). 
The internal energy equation only depends on the second order of 
the gravitational potential, through the coupling with the gas 
velocity field, ~ Vv ~ V 2( £, and therefore in this way the bulk of 
the (possible) N-body noise heating should be greatly reduced. 

We performed two re- simulations of R0 adopting M t hr = 1.1 
(R12) and M t hr = 1.5 (R13) (long dashed lines in the right panel 
ofFig|3. 

We verified the difference in the total number of cells ad- 
vanced with the total energy equation in the standard switch (772 
condition) or with the new switch {Mthr condition) by applying 
the two methods to the same snapshot of the R0 run, at z = 2 
and z = 1. We found with the first approach the number of cells 
advanced with the total energy equation is ~ 70 (~ 63) per cent 
of the total in the AMR regions at at z = 2 (z = 1); with the 
Mthr — 1-1 condition this ratio is decreased to ~ 50 (~ 45) per 
cent of cells at z — 2 (z = 1), and with M t hr = 1.5 the ratio is 
- 40 (~ 30) per cent of cells at z = 2 (z = 1). 

When the switch conditions on the Mach number are ap- 
plied since the beginning of the cluster simulation, the run with 
Mthr = 1.5 (run R13) shows a significantly reduced (by a fac- 
tor no 50 per cent) entropy production at all radii from the clus- 
ter center, but yet producing a very flat entropy distribution for 
r < 0.1R V i r . However run R13 shows a nearly constant decrease 
of entropy at all radii from the cluster center compared to all the 
other runs. Considering that the distribution of thermalized energy 
at shocks in our galaxy clusters simulations is a very steep func- 
tion of M, with a well defined peak around M ~ 2 (e.g. Vazza, 
Brunetti & Gheller 2009; Vazza et al.2010), the filtering procedure 
of run R13 is expected to filter out a ^ 50 per cent of the energy 
input from cosmological shocks. Therefore this explains why the 
large value of M t hr = 1.5 removes a significant part of the gen- 
uine production of entropy at cosmological shocks, and suggests 
that the entropy level found in the other run is mainly produced by 
the action of physical shocks dissipation. 

When the Mach number for the switch is set to the same value 
used by Springel (2010), M thr = 1.1 (run R12), the final entropy 
profile of our cluster is found to be nearly identical to the fiducial 
run, with a very similar flat entropy core. This is somewhat at vari- 
ance with the findings reported in Springel (2010), where spurious 
entropy production was masked out by applying the M t hr = 1.1 
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Figure 8. Comparison of the profiles obtained for different choices of 
the gravitational softening adopted in the N-body calculation (solid lines: 
e soft — 50kpc/h; dotted lines: e so f t = 25kpc/h) for run adopting the 
standard refinement (R2,R3) or the velocity based one (R0,R1). The vertical 
arrows show the softening length of the various runs. 



condition. A possible explanation for this may be related to the 
different on-the-fly strategies adopted to compute in run-time the 
Mach numbers between cells; further tests (also using other grid- 
based codes) may help to understand at which extent the exact for- 
mulation of the Mach number switch can be crucial to monitor the 
production of artificial entropy in cosmological cluster simulations. 

To summarize what provided by our tests with ENZO AMR, 
the choice of the flow Mach number to decide whether or not the to- 
tal energy must be used to evolve the gasdynamics of cells, is found 
to affect the normalization of the cluster entropy profile at z = 0, 
but not the shape of the distribution (at least for the M t hr ^ 1.5 
cases examined). Therefore the flat entropy core seems a pretty 
stable feature produced by ENZO AMR runs, and does not seem 
strongly affected by the details of the formalism adopted for the 
dual energy equation. At the present stage, we consider unlikely 
that the numerical issues discussed above can be fully responsibile 
for the well known difference in the entropy profiles measured ac- 
cording to SPH or to grid methods. 



3.4 The role of the smoothing in the gravitational force. 

ENZO uses a particle-mesh N-body method (PM) to follow the dy- 
namics of collision-less systems (e.g. Hockney & Eastwood 1981). 
DM particles are distributed onto a regular grid using the cloud-in- 
cell (CIC) interpolation technique, forming a spatially discretized 
DM density field. The DM density is then sampled onto the grid and 
the baryon density (calculated in the hydro method of the code) is 
added, and the gravitational potential is calculated on the periodic 
root grid using Fast Fourier Transform algorithms to solve the ellip- 
tic Poisson's equation. Since the acceleration is the gradient of the 
potential, the values of two potentials in close cells are required to 
calculate it, and thus the effective force resolution is about twice of 
the cell size. To calculate more accurate potentials on sub grids in 
the case of adaptive mesh refinement, the DM distribution is resam- 
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pled onto the finer grids using the same CIC scheme as for the root 
grid. The new boundary values are obtained with the interpolation 
from the gravitational potential on the parent grid, and a multirelax- 
ation technique is used to obtain the gravitational potential at every 
point within the sub grids (e.g. O'Shea et al.2005). 

The softening length, e so ft, used to compute the gravitational 
force and to update the DM particle motions is bound to be a mul- 
tiple of the cell size. In ENZO e so ft can be as small as the finest 
refinement level in the volume, or an integer multiple of this. 

Since the size of the entropy core is usually not much larger 
than the typical values of e so ft found in the innermost cluster re- 
gion of most of AMR simulations, it is possible that the entropy 
core is an artifact due to the fact that for scales smaller than e ao /t, 
the gravitational force does not longer obey the oc 1/r 2 scaling, but 
the oc l/(r + tsoft) 2 scaling instead0. 

We investigated this issue by producing a several runs with 
identical initial but varying the minimum allowed e so ft to compute 
gravity forces. For the the velocity-based refinement, we compared 
the fiducial run RO with run Rl (e so ft = 25kpc/h) and with run 
R19 ( e soft — lOOkpc/h); for the standard refinement scheme, we 
compared run R2 with run R3 (e so ft = 25kpc/h), with run R18 
feso/t = lOOkpc/h) and with R20 (e ao /t = 12.5kpc/h). Figure 
[8] shows that for a softening length of e so ft ^ 50kpc/h the flat 
entropy profile is a well converged feature in this cluster run, within 
a scatter of a 10 per cent at most, which may be due to slightly 
different timings in the different runs. For all softening smaller than 
tsoft = 50kpc/h, we also confirm that the two mesh refinement 
scheme always produce well separated entropy floor inside r < 
200kpc/h. For the sake of completeness, the same analysis has 
been repeated in the case of a major merger cluster (see Appendix), 
leading to consistent results. 

Mitchell et al.(2009) performed similar tests using FLASH, 
and also reported that the resolution used to compute gravity forces 
does not play an important role in setting the entropy profile inside 
the cluster, at least for maximum e so ft < 40/cpc. 

In general, we conclude that the flat entropy core in AMR sim- 
ulations is not due to spatial undersampling of the gravity forces, 
and that the flat entropy core is observed even when e so ft ~ 
0.05 — 0.1r CO re (as in run R20). We thus conclude that the en- 
tropy core in grid simulation is not dependent on the adopted grav- 
itational softening, for minimum softening length of 50kpc/h or 
smaller. 



3.5 The role of gas resolution. 

The maximum cell resolution adopted to compute the hydrodynam- 
ical equations may play a role in setting the entropy content of a 
galaxy clusters, modifying the dynamics and propagation of shocks 
waves and modifying the evolution of chaotic motions in the ICM. 
The distribution of shocks energy and Mach number in grid codes is 
a quite regular function of the underlying grid resolution (e.g. Ryu 
et al.2003; Vazza, Brunetti & Gheller 2009). The average Mach 
number responsible for most of the thermal energy dissipation is a 
well converged quantity, for grid resolutions better than ~ 500/cpc, 
even if the exact distribution of the high Mach number (M > 10) 
tails is dependent on the underlying grid resolution. 

However in the tests usually found in the literature for grid 

3 We note that however in recent ENZO AMR re- simulation of more mas- 
sive galaxy clusters the scale of the entropy core is found to be ~ 10 larger 
than the gravitational softening (Vazza et al.2010). 
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Figure 9. Overview of the effects of gas resolution and gravitational soft- 
ening in the re- simulations. Red lines: runs adopting the velocity-based re- 
finement (R1,R4,R6). Blue lines: runs adopting the overdensity-based re- 
finement (R3,R5,R7). The additional black lines are for the fiducial run (R0, 
solid line) and for run R2 (dotted line). The vertical arrows show the soft- 
ening length of the various runs. 



simulations the above dependence is coupled with the dependence 
on e so ft (which may also indirectly affect the production of shocks 
by changing the morphology of in-falling gas/DM sub-clumps). We 
aim here to disentangle the two effects, by keeping the maximum 
softening length fixed to e so ft = 25kpc/h, bud varying the maxi- 
mum gas resolution in the computation of fluid-dynamical effects. 

In runs Rl and R3 we allow ENZO to refine up to one level 
more, I = 4, keeping the same setup and refinement strategy of 
runs R4 and R5, respectively. If we compare the solid lines in 
Fig|9](run R1,R3) and the dotted lines (run R4,R5) we see that 
the trend with gas resolution is opposite in the two mesh refine- 
ment strategies: while in the velocity-based refinement the entropy 
inside r ~ lOOkpc/h is increased by a ~ 20 — 30 per cent (R4), 
in the case of the overdensity-based refinement the inner entropy 
is decreased by a similar amount. The gap in entropy inside r cor e 
is of the order of AS ^ 30 — 60keVcm 2 . We also re-simulated 
runs R4 and R5, allowing the code to increase also the softening 
length up to I = 4 (run R6, R7, long dashed lines): the reported 
trends are the same (see red and blue dashed lines in Fig(9}. The 
re- simulations shows that the difference between the two refine- 
ment schemes is always significant: a very flat entropy core inside 
r < lOOkpc/h is produced in both cases, but a significant gap is 
found when comparing the two strategies, with the implemented re- 
finement scheme producing the larger value. The trend is confirmed 
also by the comparison of two runs where the minimum softening 
and cell resolution were fixed to I = 2 (50kpc/h), run R22 and run 
R21 (dot-dashed lines in the same panel). 

The reason for the opposite trend in the velocity-based scheme 
can be understood by the increased role of mixing motions, which 
cause a slightly more efficient inward transport of higher entropy 
material from the cluster outskirt (see also SecHJ). The efficiency 
of this mechanism is expected to increase as resolution is in- 
creased, since the velocity-based strategy is explicitly designed 
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Figure 10. Radial entropy distribution for all TV = 10 5 tracers evolved 
in run RO. The colored sub sample shows the selections used for tests in 
Sec |4.1| The additional lines show the average entropy profile of the clus- 
ter (red/dashed:gas density weighted profile; orange/dot-dashed: volume- 
weighted profile). 



to reduce the artificial dampening of small scale chaotic motions, 
which are partly responsible of diffusive mixing in the ICM. In the 
overdensity-based approach clumps are refined more and more, and 
they can deliver low entropy gas in the innermost region in a more 
efficient way; an increase of resolution also minimizes the effect of 
numerical mixing and let the cold gas phase to survive longer (e.g. 
Wadsley et al.2008; ZuHone, Markevitch & Johnson 2009). On the 
other hand in the velocity-based strategy the clumps are also more 
refined, but they are also more efficiently destroyed, before reach- 
ing the cluster center, by ram-pressure stripping at the outer regions 
and excite chaotic motions which are more long-living since they 
are not damped by the code; the net effect is an increase of entropy 
inside r cor e in the velocity-based strategy. This stresses the need of 
having an accurate description of shocks and turbulent motions of 
the ICM, since the inner entropy distribution does not only depend 
on the maximum resolution within the core, but also on the reso- 
lution at the outer regions, where the bulk of the cluster entropy is 
produced. Convergence tests reported in Mitchell et al.(2009) and 
in ZuHone, Markevitch and Johnson (2009) suggest that full con- 
vergence for the entropy profile in the standard refinement strategy 
is reached for a maximum resolution of ~ lOkpc/h or smaller. 
A similar conclusion is likely also the velocity-based strategy, but 
we could not run so far a re- simulations reaching / = 5 for com- 
putational limitations. However the trend of the red lines in Fig|9] 
suggest that convergence is near. 

We conclude that even if the different mesh refinement strate- 
gies and the adopted maximum refinement levels can produce mod- 
ification in the level of the inner entropy budget (up to a factor of 2 
at the peak resolution investigated), the flat entropy core is present 
in all cluster resimulations (see also the additional resolution tests 
in the Appendix). This again calls for a mechanism of physical na- 
ture, which drives an efficient spreading of entropy enriched gas in 
the innermost regions of evolved clusters; this will be investigated 
in detail in the next Section. 






redshift 

Figure 11. Time evolution of gas density (top panel), gas entropy (central 
panel) and mean gas velocity (or mean gas turbulent velocity, dashed lines) 
for the 5 groups of tracers selected according to their position in the cluster 
atmosphere at z = 0, as in Fig[l0l(see Sec|4]for details). 
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Figure 12. Time evolution of gas density (top panel), gas entropy (central 
panel) and mean gas velocity (or mean gas turbulent velocity, dashed lines) 
for the 5 groups of tracers selected at z = 2 according to their initial gas 
density (see Sec|4]for details). 




0.01 0.10 1.00 

r[Mpc/h] 

Figure 13. Radial profile of the mixing parameters for the tracers in run R0 
(dashed lines) and radial profile of gas entropy (solid lines) , at z = 1 (blue), 
z = 0.5 (green) and z = (red). The long-dashed line shows the radial 
profile of 0.5 — M at z = 0. The gas entropy has been dived by 2000 for 
a better visualization. The vertical arrow show the softening length adopted 
in this run. 

4 PHYSICAL EFFECTS ON THE ENTROPY 
PRODUCTION. 

The numerical tests reported in the previous Sections have shown 
that the presence of a regular low gas entropy distribution in galaxy 
clusters simulated with grid-based techniques is a not a numerical 
artifact, rather but a very stable feature against a number of impor- 
tant changes in the possible setup of a cosmological simulation at 
high resolution. 

At this point it is interesting to answer to the following ques- 
tions: a) what is the main physical mechanism which sets the in- 
ner gas entropy distribution in a forming galaxy cluster, in non- 
radiative runs? b) Is the inner gas entropy distribution affected by 
a more sophisticated physical modeling of cluster dynamics (e.g. 
employing radiative cooling)? c) What is the effect of other non- 
gravitational extra-heating mechanisms (e.g. AGN feedback) on the 
the gas entropy profile of clusters? 

In the following Section we explore how entropy is build over 
time in the same galaxy cluster analyzed above, by means of a La- 
grangian approach based on tracers particles. In the other Sections, 
we analyze how is the entropy floor modified when a more realistic 
modeling of cluster physics is considered (e.g. assuming radiative 
cooling) and if extra-heating processes (e.g. uniform pre-heating or 
AGN feedback) are capable to reproduce an entropy distribution 
similar to that of the fiducial run. 



4.1 Where and when entropy is build in a forming cluster. 

Passive tracers are a useful tool to follow the average trajectory 
of accreted baryons in a growing galaxy clusters, and to study the 
mixing pattern driven by accretion phenomena in clusters (Vazza, 
Gheller & Brunetti 2010). Here we use tracer particles to track the 
exact origin in time and space of the gas entropy deposited in the 
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Figure 14. Comparison of the entropy profiles of radiative runs, using the 
standard mesh refinement (blue) and the velocity-based refinement (red). 
The additional dashed lines show the profile of the corresponding non- 
radiative runs (R0 and R2). The vertical arrow shows the softening length 
adopted in run R15 and R16. 



cluster core of run R0, and to explain the emergence of a flat en- 
tropy core. We inject passive (mass-less) tracers in the simulation 
and let them be passively advected in the ICM using the informa- 
tion of the 3-D velocity stored for each snapshot of the simulation, 
in a post-processing phase. Convergence tests on the interpolation 
technique to assign velocities to tracers, on the time and space sam- 
pling are discussed in detail in Vazza, Gheller & Brunetti (2010). 

In this particular case, we injected N = 10 5 tracers at z = 2 
in the fiducial run (R0), within a cubic volume of « (10Mpc/h) 3 
centered in the AMR region. At this time the total virial mass of the 
forming cluster is only 10 — 15 per cent of its final mass at z = 
(Fig[TJ and only small proto-clusters are present within the AMR 
region. 

Tracers are initially placed with a random sampling of the vol- 
ume, and then their positions are updated with the time steps finely 
saved in time of the original simulation, using a Nearest Grid Cell 
interpolation scheme (e.g. Hockney & Eastwood 1981). At all time 
steps, the tracers record the thermodynamic values of the closest 
cell in the grid distribution, and the whole thermodynamic history 
along the trajectory of every tracers can be recovered for analysis. 

In Figure [TO] we show the radial distribution of tracers at 
z = 0, plotting on the vertical axis the entropy of the nearest cell at 
each location. As a comparison, we overplot the (density weighted 
and volume weighted) entropy profile of the cluster, to confirm 
that tracers sample the underlying Eulerian distribution in an ac- 
curate way. At z = we selected 5 shells of tracers with width 
25kpc/h, located at the radii of r < 50kpc/h, r = 150kpc/h, 
r = 300kpc/h, r = 600kpc/h and r = 1200kpc/h (shown as 
different colors in FigfTUt. This allowed us to study the mean evo- 
lutive history of the parcels of gas ending at the different level in 
the cluster entropy profile at the end of the simulation. In FigfTTlwe 
show the behavior of the mean gas entropy, gas entropy and veloc- 
ity (or chaotic velocity, as measured in Sec l3.2l ) modulus for the 5 
different shells, as a function of evolving cosmic epoch. Except for 



the first bin (which contains N — 20 tracers) for all the other shells 
the number is of the order of ~ 200 and thus the mean values are 
very robust. 

The sharp spikes in gas density (top panel) trace the main 
merging episodes which involve the different "shells" of tracers. 
The cluster entropy profile at z = is manifestly produced by a 
uniform and regular in time mechanism of "sorting in entropy", 
which affects every tracers at the moment of its entrance in the 
virial region of the forming cluster (middle panel). The level of gas 
entropy in the cluster core is set, on average, during intense shock 
heating at z ~ 1 (which also corresponds to the epoch of the most 
net increase in cluster mass, as shown in Fig [I). On the other hand 
the gas tracers with a larger entropy (S ~ lOOOkeV cm 2 ) at the 
final epoch are found to be intensely shock heated at more recent 
epochs, z rsj 0.6 — 0.8. This analysis show that on average they 
belong to smooth, low density gas environment at z = 2, and that 
they are subject to smaller bulk and chaotic velocity fields at the 
final epoch, since they do not belong to bound in-falling structures 
(lower panel). 

A complementary test was run by directly selecting 5 families 
of tracers at z = 2 and sorting them according to their initial gas 
overdensity; their evolution is followed in Fig[l2] The "sorting in 
entropy" among tracers is even more evident with this setup: the 
average values of gas density and gas entropy for the different fam- 
ilies never overlap for z < 1, meaning that the main variable which 
sets the final entropy of a gas parcel ending up into a massive galaxy 
cluster is its initial overdensity. Say it differently, the fact that a gas 
particle is in a overdense (clumpy) environment around the form- 
ing cluster, determines on average the time at which it gains the 
bulk of the final entropy (in its first impact on the volume of the 
cluster under virialization), and its final distance from the cluster 
center, through the mechanism of entropy sorting in the main clus- 
ter atmosphere. This is in excellent agreement with the spherical 
analytic models of clusters forming in a hierarchical scenario, that 
prescribes a raising entropy at increasing radius from the cluster 
center, following the progressive deposition of shells undergoing 
stronger and stronger shock heating during the hierarchical growth 
of a cluster, (e.g. Tozzi & Norman 2001; Cavaliere, Lapi & Fusco- 
Femiano 2009). 

This simple scenario may of course be more complicated in 
the case of a major merger, which can mix the intra cluster medium 
in a more efficient way (e.g. Mc Carthy et al.2007; Poole et al.2008; 
Vazza, Gheller & Brunetti 2010; ZuHone 2010). 

The above tests show where and when the entropy profile is 
produced during the cluster evolution, but do not necessary imply 
the emergence of a small inner region of size r cor e ~ 0.1R V i r 
where the stratification is broken and the gas sets to the constant 
value of S ~ lOOkeVcm 2 observed in the previous runs. A viable 
mechanism naturally produced by the accretion of matter onto the 
cluster is mixing of the inner gas layers, in response to chaotic mo- 
tions in the ICM. To better compute the degree of mixing between 
Lagrangian tracers after their injection at z = 2, using a formal- 
ism introduced in Vazza, Gheller & Brunetti (2010). We introduce 
the mixing parameter of an "s" family of tracers, M s , respect to all 
existing "i" families, by computing: 

\n s — V\ m(i ^ s)\ 
Ms = l- 1 — % n t. >\ (2 ) 

where n s is the number density of the "s" tracers within a cell 
(at the highest resolution level) and the sum refers to all the species 
of "i" tracers (included the "s" specie). This formula generalizes 
the more simple case of mixing between two species (e.g. Ritchie 
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& Thomas 2002) and has a simple interpretation: for a cell where 
m « n<2 « 77,3... w n s the different families are well mixed and 
we have M s — >• 1, while M s — >• implies no mixing within the 
cell. 

The total mixing in each cell, M, is the volume average be- 
tween all species, Mj = M S /N SJ where iV s is the number of 
families considered. To obtain better statistics, the cell considered 
to compute Mi have a size which is double than the maximum gas 
resolution (50kpc/h). 

We sorted the density of all tracers at z = 2 and divided the 
sample in 5 bins, each with 2000 tracers; then we computed the 
number density of each "s" family of tracers within the cluster, and 
computed the mean mixing Mi at different redshifts (Fig [13]), only 
for those cells containing at least 1 tracer. 

The profile of the mean mixing presents a regular trend in 
time, with a inner core of mixing material which is build as the 
cluster accretes matter, until an almost flat profile (with mean mix- 
ing M % & 0.3) is found for r < lOOkpc/h at z = 0. Outside of this 
radius, Mi has a sharp decline towards the virial radius of the clus- 
ter, where the mixing is very poor Mi < 0.05. We verified that also 
different initial choices for the setup of the tracers sampling (e.g. 
by adopting a different number of tracers, or a different number of 
species) do not affect the above trend in any significant way. 

FigEl shows the evolution of the gas entropy profile (solid 
lines) and of the mean mixing profile (dashed lines) for three 
redsfhits. Quite clearly, the formation of the central entropy floor 
happens together with the formation of the prominent mixing pat- 
tern in the center of the cluster (the similarity of the two profiles 
is better indicated by overplotting the radial profile of | — Mi | at 
z = 0, long dashed line), and the S(r) oc r 1 " 1-1 " 2 scaling ex- 
pected from SPH simulations (Voit et al.2005, Fig[6l is broken 
starting from the same radius for which Mi sharply increases to 
its maximum in the center. In Vazza, Gheller & Brunetti (2010) 
we provided evidence that in general the building of cluster en- 
tropy and of prominent mixing patterns are correlated features 
also in major merger systems (where, however, large and transient 
plumes of efficient mixing can be found also at large cluster radii); 
also recent FLASH AMR simulations (Zu Hone 2010) has lead to 
similar conclusions, based on a parametric study of binary clus- 
ter mergers. Usually the innermost region of cosmological cluster 
runs is characterized by small scale subsonic motions (on scales 
< 500kpc/h and a v ~ 0.3 — 0.5c s , where c s is the gas sound 
speed) continuously excited by the crossing of gas/DM material 
accreted within the cluster (e.g. Norman & Bryan 1999; Dolag et 
al.2005; Iapichino & Niemeyer 2008; Vazza et al.2009). In the case 
of ENZO AMR simulations these motions are well characterized 
by a power law spectrum for nearly two orders of magnitude in 
spatial scales (Vazza et al.2009; Xu et al.2009; Vazza, Gheller & 
Brunetti 2010); also the tracers pair-dispersion statistics show a 
well defined power-law dependence on time (P(t) ~ t 3//2 , where 
P(t) is the distance between couple of tracers initially located at 
a small distance) compatible with a fairly fast transport motions in 
the turbulent ICM (Vazza, Gheller & Brunetti 2010). 

When combined together, all the above evidences lead to the 
consistent conclusion that the emergence of a flat entropy floor in 
clusters simulated with grid codes is mainly due to the integrated 
effect of (mostly subsonic) turbulent motions in the evolving ICM. 
These motions are quite effective in mixing parcels of gas, which 
have been initially sorted according to the "entropy sorting" mech- 
anism described above, generally inside r < 0.1R V i r and for late 
z < 1 redshift, when that the bulk of cluster mass has been assem- 
bled. The same mechanism also applies when the standard mesh re- 



finement is adopted, provided that the turbulent energy in the clus- 
ter core is smaller (Fig[5j and that the mixing and the transport of 
gas particles is significantly smaller (e.g. Fig. 17 in Vazza, Gheller 
& Brunetti 2010). 

It is well known that grid codes are prone to numerical mix- 
ing (i.e. different gas phases are forced to combine into an average 
cell value when their separation is smaller then the minimum avail- 
able cell size) while in SPH gas particles do not mix entropy by 
construction, unless ad-hoc diffusion term is considered in the SPH 
equation (e.g. Agertz et al.2007; Price 2008; Wadsley et al.2008; 
Merlin et al.2010). Numerical mixing must be considered as an 
additional source of mixing also fro the ENZO AMR resimula- 
tions presented here. However the ubiquitous finding of evident 
mixing motions on scales much larger than the cell size discussed 
above, combined with the evidence that the entropy floor presents 
a very small evolution with the grid resolution, for peak resolution 
< 25kpc/h, indicate that numerical mixing cannot be the main 
responsible for the production of the entropy floor. According to 
this interpretation, the entropy floor is thus a result of mixing of 
the baseline steep S(r) ~ r 1A profile (produced by the progres- 
sive shock heating of infalling shells of matter) within the region 
where the turbulent energy is maximum within cluster, < 0.1R V i r , 
as a result of the continuous excitement of a hierarchy of chaotic 
motions in the ICM driven accretion after accretion. 

This is true for an unviscid treatment of the ICM, while the 
presence of magnetic fields and plasma viscosity may alter this pic- 
ture in a significant way (e.g. Parrish & Quataert 2008; Ruszkowski 
& Oh 2010 and references therein). 

In the following Sections, we will explore more realistic mod- 
elization of the same cluster, where non-gravitational mechanisms 
of entropy decrease, such as radiative cooling, and entropy in- 
crease, such as heating from energy feedback by stars of active 
galactic nuclei, are computed in run time. 

4.2 The role of radiative cooling and pre-heating 
mechanisms. 

Radiative cooling in galaxy clusters have a dramatic impact in the 
thermodynamics of the ICM, if no heating mechanism other than 
shock heating is available to balance the cooling catastrophe, for 
those cluster regions in which the gas cooling time is < < of the 
cluster age (e.g. Katz & White 1993). Figure [T4l illustrates the ra- 
dial profiles of gas density, gas temperature and gas entropy for the 
two re- simulations with radiative cooling (assuming a fully ion- 
ized H-He plasma with a constant metallicity of Z — 0.3Z©) and 
with standard mesh refinement (R15) or velocity-based refinement 
(R16). In both cases, a steeply decreasing entropy profile devel- 
ops towards the center of the cluster for r < lOOkpc/h, with very 
low entropy values, S « IkeVcm 2 in the center, a massive gas 
condensation peaking at 10 -24 gr / cm 3 and a temperature dip with 
T « O.lkeV, similar to the classic theoretical cooling flow sce- 
nario (e.g. Fabian 1994). 

We note that these simulations do not consider any prescrip- 
tion for star formation from the dense and cold phase of the ICM, 
and therefore the central condensation produced by run R15, R16 
(see the solid black line in the top panel of FigfT6t would be sig- 
nificantly reduced by modeling star formation in a self-consistent 
way (e.g. Pearce et al.2000; Yoshikawa, Jing & Suto 2000; Val- 
darnini 2002). However this would make the pure cooling runs and 
the other explored in the next Sections more computationally ex- 
pensive. We thus preferred to defer to the future the study of this 
issue, and to use this simplified cooling model without star forma- 
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Figure 15. Maps of gas density (top panel) and gas entropy (bottom panel) 
for a slice taken in run B2 (cooling, pre-heating and refinement on velocity 
jumps). The side of the image and the meaning of colors are as in Fig 0) 



tion as the framework to study the effects of thermal energy feed- 
backs in the ICM, and the modifiation they cause to the gas entropy 
distribution. 

In general major and moderate mergers (e.g. Burns et al.2008; 
Poole et al.2008) or gas sloshing triggered by the passage of 
DM/gas sub-clumps (e.g. Ascasibar & Markevitch 2006; ZuHone 
et al.2009) may significantly reduce the cooling catastrophe in ra- 
diative simulations, by exciting internal merger shocks or turbulent 
mixing in the clusters core. However, in the particular cluster simu- 
lated here the amount of chaotic motions excited for z < 1 is never 
powerful enough to slow down the cooling flow in any significant 
way, even when the additional mesh refinement criterion is turned 
on. We report that also in the case of the merger system studied in 
the Appendix, the major merger at z ~ 0.85 does not reduce the 
cooling catastrophe at the end of the simulation. 
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Figure 16. Gas density and entropy radial profiles for the run adopting pre- 
heating (run PHI So = lOkeVcm 2 , run PH2 So = lOOkeVcm 2 and run 
PH4 So = 200/ceVcra 2 ). The dotted lines show the result for the run PH3, 
adopting So = lOOkeVcm 2 and standard mesh refinement. The additional 
black lines shows the profile for the non-radiative run (run R0, long dashed) 
and for the cooling run (run R15, solid). We also overplot in blue the best fit 
profiles of Cavagnolo et al.(2009), with a core entropy of So = lbkeVcm 2 
(dot-dashed) and So = 150keVcm 2 (long-dashed). The vertical arrow in 
both panel show the softening lenght adopted in the runs with cooling and 
pre-heating. 



As widely known, cluster configurations according the "pure 
cooling" scenario are not observed in the real Universe, and the 
gas temperature in real clusters is never observed below ~ O.lkeV 
(e.g. Rossetti & Molendi 2010 and references therein). For this rea- 
son, additional sources of gas heating were considered in order to 
reconcile simulations with observations. In what follows, we will 
apply some of the most promising models of extra-heating devel- 
oped in the literature to our cluster simulations, and study their im- 
pact on the gas entropy profile at z = 0. 
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4. 2. 1 Early pre-heating. 

Motivated by the early evidences of significant departures from 
self-similar scalings expectations in observed galaxy clusters (e.g. 
White 1991; Evrard & Henry 1991; Ponman , Cannon & Navarro 
1999; Lloyd-Davies, Ponman & Cannon 2000), several authors 
proposed a "pre-heating" scenario, in which an energy input of 
non-gravitational origin is injected in the Intra Galactic Medium 
at early redshift (z ~ 3 — 10). This smoothens the gas of in-falling 
primordial halos of low mass and supplying the innermost region 
of massive cluster of an amount of high entropy gas, to recon- 
cile with observations, (e.g. White 1991; David, Forman & Jones 
1991; Kaiser 1991; Evrard & Henry 1991; Cavaliere et al.1997; 
Voit et al.2005). This early heating may be provided by a num- 
ber of sources, such as star formation and SNe explosion, radiative 
and mechanical heating from AGNs, radiative heating from hard 
X-ray background, etc (e.g. Tozzi & Norman 2001). The estimated 
needed amount of extra entropy at high redshifts falls in the range 
100 ^ AS ^ 300keVcm 2 , and viable sources of it in the early 
Universe can be supernovae explosions, star formation and galactic 
winds. From the numerical viewpoint, several group attempted to 
model this process either by imposing an impulsive entropy injec- 
tion at a given epoch (e.g. Bialek, Evrard & Mohr 2001; Borgani et 
al.2005; Kay et al.2007; Romeo et al.2006; Younger & Bryan 2007) 
or in a redshift-modulated way (e.g. Borgani et al.2002; Sijacki et 
a.2007; Mc Carthy et al.2009). 

An important drawback of early pre-heating models could be 
that they tend to remove to much low entropy gas from lower mass 
halos, without obtaining a realistic galaxy population (e.g. Don- 
ahue et al.2006). In addition, recent results based on XMM-Newton 
analysis presented by Rossetti & Molendi (2010) proved that most 
of non-cool core clusters host regions with low entropy and high 
metallicity, suggesting the possibility of a recent transition between 
cool core and non-cool core systems, contrary to the pre-heating 
scenario (e.g. Guo & Mathews 2010). 

As a first step to investigate the role played by non- 
gravitational heating on the entropy level of our cluster, we tested 
early pre-heating models following the prescription of Younger & 
Bryan (2007). 

The thermal energy within each gas cell in our cluster run is 
selectively increased at z = 10 by: 

k B AT = So[^r-\ (3) 

In detail we re- simulated run R15 (velocity based refinement 
and cooling) by imposing the entropy floor of So = lOkeVcm 2 
(run PHI), So = lOOkeVcm 2 (run PH2) and S = 200keVcm 2 
(run PH4). Smoothing the gas density distribution of halos at high 
redshift is expected to affect the shock heating process in form- 
ing structures, through "entropy amplification" at strong acrretion 
shocks (Voit et al.2005). To highlight the dependence of this mech- 
anism on the accuracy with which accretion shocks are modeled in 
our runs, we also tested the intermediate pre-heating scenario of run 
PH2 (So = lOOkeVcm 2 ) in the standard mesh refinement strategy 
alone (run PH3). 

Figure Q3] shows the effect of pre-heating in extreme case of 
run PH4 at z = 0. Compared to the non pre-heated scenario (run 
R0) the global cluster morphology is smoother, and the accretion 
pattern are more regular in shape, since most of the accreted gas 
clumps were smoothed by the early heating episode. In FigfTolwe 
compare the gas density and gas entropy profiles for the runs em- 
ploying pre-heating, against the profiles of the fiducial run (R0) and 
of the simple cooling run (run R12). 



Consistently with the literature, we find that energy inputs cor- 
responding to So = 100 — 200keVcm 2 are capable of keeping the 
in-falling gas on an higher adiabat, preventing the core gas to cool 
below ~ 0.5keV. In particular, the pre-heating prescription of run 
PH4 is suitable to recover the gas density and the gas entropy of 
the non-radiative run (R0), within a ~ 10 — 20 per cent at all radii. 
A similar result is also found for the major merger system stud- 
ied in the Appendix. On the other hand, the uniform pre-heating 
model with So = lOkeVcm 2 is found to be insufficient to prevent 
the cooling catastrophe and the resulting cluster profile is almost 
identical to that of run R15. 

The trend found is qualitatively in agreement with the re- 
sults of Borgani et al.(2005) and Younger & Bryan (2007), even 
if our study is based on a single object and no conclusion about 
the most suitable values of So needed to reconcile with observa- 
tions can be derived in a statistical sense. We also note that the 
final entropy configurations of run PH2, PH3 and PH4 are compat- 
ible with the bimodal distributions of entropy profiles of obtained 
with CHANDRA (Cavagnolo et al.2009, additional blue lines in 
FigfToT). This work recently suggested the existence of two broad 
population of clusters, characterized by an inner entropy value of 

5 ~ 15keVcm 2 or S ~ 100 — 150keVcm 2 , and a large radial be- 
havior scaling as S(r) oc r 1 ' 1 . Our cluster is a relaxed one at z = 0, 
and the fact that it is more similar to the "low entropy core" class 
of CHANDRA clusters, even when So = 200keVcm 2 is applied, 
is fully compatible with the idea that "high entropy core" class is 
produced only by those clusters with a sufficiently violent merger 
in their past (e.g. Rossetti & Molendi 2010). 

4. 2. 2 Heating from AGN jets. 

The most successful models to achieve a balance with radiative 
cooling during the simulated evolution of clusters rely on heating 
by outflows from an AGN hosted by the central massive galaxy 
(e.g. Churazov et al.2000; Bruggen & Kaiser 2002; Brighenti & 
Mathews 2003; Dalla Vecchia et al.2004; Zanni et al.2005; Heinz 
et al.2006). 

The self-consistent modeling of AGN heating in evolving 
galaxy clusters, in connection with the matter accretion history 
of the central massive cD galaxy has become only quite recently 
within the capability of full cosmological simulations (e.g. Si- 
jacki & Springel 2006; Dubois et al.2010). In general, this issue 
is made complex by the (still unclear) role played by other phys- 
ical mechanisms relevant to the thermodynamic evolution of the 
cluster plasma: plasma viscosity, magnetic fields, Cosmic Rays, etc 
(e.g. Bruggen & Kaiser 2001; De Young 2003; Brighenti & Math- 
ews 2003; Ruskowski et al.2007; Sijacki et al.2008; Xu et al.2008; 
O'Neill & Jones 2010; Scannapieco & Bruggen 2010). 

In this Section, we study the energy budget necessary to 
quench the cooling catastrophe developed in our radiative run 
(R15), using the thermal feedback from an assumed central AGN 
(e.g. Booth & Schaye 2009; Teyssier et al.2010) This is done by in- 
jecting localized inputs of extra- thermal energy, starting at z = 1 in 
run Rl 5, in the region of the cluster density peak, where the AGN 
is assumed. 

This is motivated to mimic the thermal feedback response 
from a central AGN which releases as feedback a fraction of the 
rest mass energy of the accreted cold gas, with an efficiency which 
is generally assumed to be in the range of ~ 0.01 — 0.1 of the 
total energy radiated by the super massive black hole (e.g. Booth 

6 Schaye 2009; Giodini et al.2010 and references therein). In our 
setup, we assume that the feedback manifests itself at the scale of 
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Figure 17. Cuts through the center of the simulated cluster, showing the evolution of the bolometric x-ray luminosity (top panels, arbitrary units) and gas 
entropy (bottom panels) for run B2. From the left to right, the epoch shown are z = 0.96, z = 0.85 and z = 0.75. The side of the images is 2.7Mpc/h. 



the maximum gas resolution in the cluster center (25kpc/h), as the 
injection of two point-like over-pressurized "bubbles", produced by 
the interaction between the launched jets and the surrounding cold 
ICM. The thermal energy at the location of the "bubbles" is updated 
as: 



3k B pTAx 3 
2 /imp 



+ 



(4) 



where fi is the mean molecular mass, m p is the proton mass, 
where Ax 3 is the volume of the cell, and e jet is the fraction of 
ejet released at every injection episode. The 3-D velocity field at 
the injected bubbles is left unchanged. Therefore in our treatment 
€j et represents the thermal energy released in the ICM by the the 
two jets at a given time step, after the thermalisation of a part of 
their mechanical energy, which is assumed to happen on a sub- 
grid scale. The "bubbles" are initially located at the distance of 
djet = 50kpc/h, at two opposite sides of the gas density peak 
of the cluster, starting from z — 1 in run R15 (or run R16 for the 
standard refinement strategy). The injection is performed before the 



hydro step of the PPM scheme, by updating the gas internal energy 
following EqHlof two cells (at the maximum resolution level) ac- 
cording to [4j then the Riemann solver in ENZO is evolved in the 
usual way. 

Preliminary tests showed that the impulsive injection of the 
whole ejet in a single time step of the simulation produces un- 
realistically strong (e.g. M > 10) shock waves in the cold and 
dense cooling flow cluster region. Since only mild shocks are ob- 
served in jets/bubbles interaction with the ICM of real clusters (e.g. 
Simionescu et al.2009; Werner et al.2010), we preferred to adopt 
a more gradual release of energy from the central gas peak, by 
distributing 6j e t in ~ 20 injection episodes, across a total time 
of ~ SGyr, preserving the same orientation for the "jets" axis. 
In principle, idealized but more self-consistent recipes to link the 
feedback energy with the matter accretion rate within the cooling 
region can be applyed to cosmological simulations (e.g. Sijacki & 
Springel 2007; Booth & Schaye 2009); however here we want to 
investigate how different re- simulations of the same object react 
to a constant model of extra thermal energy release from a cen- 
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Figure 18. Zoomed view of the internal ICM region after the first AGN 
outflows in run B2 at z = 0.96 (logarithmic view of gas entropy). The 
isocontours are drawn only for gas with S > 500keVcm 2 to highlight 
the contribution from the AGN burst. The scale of the images is ~ 1.5 x 
O.bMpc/h. 

tral AGN. We defer to the future any study of more self-consistent 
setup of the feedback energy, and of the way it is released within 
the ICM (e.g. by varying the jets orientations in time, by assuming 
an "quasar mode" and "radio mode" feedback, etc). 

With our setup, only M < 5 shocks are produced, even in the 
most extreme scenario investigated, and only in the starting phase 
of the jet, when the surrounding ICM is in its coldest phase. Our 
trials adopted 6j e t = 10 58 erg (run Bl) and €j e t = 10 59 erg (run 
B2); this makes the typical power of our jets in injection phase of 
about Wjet ~ e jet /tste P ~ 10 42 -10 43 erg/s (t s te P is ~ 4- 10 7 yr 
at that epoch). We note that the power for the energy release of our 
jets in the surrounding medium, and the assumed duty cycle and 
duration are within most of the estimated energy budget provided 
by the observations of AGN activity reported by many authors (e.g. 
Birzan et al.2004; Dunn & Fabian 2006; Wise et al.2007; Giacin- 
tucci et al.2008; Bird, Martini & Kaiser 2008; Worrall et al.2009; 
Liuzzo et al.2009; Sanders & Fabian 2009; Gu, Cao & Jiang 2009; 
Gitti et al.2010; Giodini et al.2010). 

Scannapieco & Bruggen (2008) have recently shown that a 
proper treatment of turbulence on < lOkpc scales is mandatory to 
model the full interaction between jet-inflated bubbles and the ICM, 
because this may change the rate of energy transfer to the surround- 
ing cold phase of the ICM. Therefore, it is unlikely that our simula- 
tions are fully converged, and further tests at higher resolution will 
be needed in the cosmological framework. In any case, we assess 
the role played by numerical resolution here by running two ad- 
ditional re- simulations, using only standard mesh refinement, (run 
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Figure 19. Gas density and entropy radial profiles for the run adopting AGN 
jets feedback (run B1,B4 Cj et = 10 58 ergs, runs B2 ej et = 10 59 ergs). 
The dashed lines show the result for the run B3, adopting Cj et = 10 59 ergs 
and standard mesh refinement. The additional black lines shows the profile 
for the non-radiative run (run R0, long dashed) and for the cooling run (run 
R15, solid). We also overplot in blue the best fit profiles of Cavagnolo et 
al.(2009), with a core entropy of So = 15keVcm 2 (dot-dashed) and So = 
150keVcm 2 (long-dashed). The additional arrow in both panel shows the 
softening length adopted in the runs with cooling and AGN feedback. 

B3) and using an additional level of refinement in the velocity based 
strategy, up to a maximum resolution of Ax = 12. 5kpc/h (run 
B4). 

The evolution of bolometric X-ray luminosity and gas entropy 
in a slice crossing the cluster center for run B2 is shown in FigfTTl 
Soon after the first injection, two vertical outflows has developed 
for ~ 400 — 500/cpc along the axis of "bubbles" injection and has 
pushed the dense and low entropy material out to larger radii. The 
first feedback episode drives a mild shock in the cold ICM, with 
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M rsj 3.5, while along the outflows Kelvin-Helmoltz instabilities 
develop and favors the mixing between the cold uplifted gas from 
the core and the surrounding hotter ICM (Fig j4.2.iT ). The outflows 
that follow inflate more stable "bubbles" (central panel in FigfTTl, 
which are initially less overpressurized compared to the surround- 
ing ICM, heated by the previous feedback episode. These bubbles 
only drive weak M ~ 2 shocks around the central gas condensa- 
tion. In addition, the sloshing motions in the cluster center are pow- 
erful enough to bend the initial orientation for the bubble launching, 
and to partially provide azimuthal mixing of the injected entropy. 

Compared to the fiducial run (RO) or to the pre-heated runs 
(PH2,PH4), we found that the large scale accretion patterns are 
modified by the outgoing propagating shocks that follow the AGN 
activity. 

In Figure [19] we compare the the gas density and gas entropy 
profiles for all trials at z — 0. We found impossible to recover an 
entropy profile similar of the fiducial run (and also with the results 
of Cavagnolo et al.2009, overplotted in the same Figure): the extra 
energy of shocks is very efficiently delivered to larger cluster radii 
by the shocks, which develop in the interaction between the cooling 
ICM and the hot bubbles phase. As a result, in these configurations 
we find an excess of gas entropy for r > 75kpc/h, compared to the 
radiative and non-radiative cases and a flat entropy profile inside the 
cluster core. However these runs produce a much smaller cooling 
region (r CO oi < 50kpc/h) compared to the r coo i ~ lOOkpc/h of 
pure cooling models (run R15, R16). The observed steepening of 
the internal gas entropy after the AGN energy released is in good 
agreement with the semi- analytical prediction of a quasar-driven 
blasts in clusters discussed in Lapi, Cavaliere & Menci (2005), fur- 
ther suggesting that in the trials investigated here the shock heating 
mechanism is the main mechanism which interchanges energy be- 
tween the AGN and the surrounding ICM. 

When the mesh refinement is triggered uniquely by gas/DM 
over-density (run B3), the final entropy profile at z = present a 
larger cooling region (in-between the pure cooling case and the B2 
re- simulation). In this case it is difficult to disentangle the effect 
of the under- sampling issues of satellites-driven mixing (Sec l4.ll ) 
from that of the under- sampling of jet-driven turbulence around the 
cluster core. The re- simulation with an additional level of mesh 
refinement (run B4) shows that as resolution is increased the final 
size of the cooling region is reduced, and the entropy at larger radii 
is increased, due to a better modeling of shocks induced by the 
outflows and of the driven turbulent motions. Further studies will 
be needed in the future to fix the best resolution needed for full 
numerical convergence in these features. 

Our general conclusion is that, even if the action of feedback 
from outflows in our simplified implementation efficiently reduces 
the size of the cooling flow region compared to a pure radiative 
run, it remains difficult to reproduce a flat inner entropy profile as 
in the non-radiative case. The problems is not in the energy bud- 
get assumed in the outflow (which is a reasonable energy budget 
available to observed AGNs) but rather in the mechanism which 
transfers to the surrounding medium, which is mainly shock heat- 
ing of the cold central ICM (see also Zanni et al.2004; Lapi, Menci 
& Cavaliere 2005). It is likely that more gentle mechanisms of 
feedback from the central AGN, such as a more gradual deposition 
of many by inflated bubbles(e.g. Churazov et al.2001; Bruggen et 
al.2007), can be more efficient in stopping the catastrophic cooling 
in our run. However in that case a significantly larger resolution 
than the one available here must be considered, which is presently 
difficult for AMR runs with our mesh refinement scheme. In addi- 
tion, sub-grid modeling of turbulence may necessary (e.g. Scan- 



napieco & Bruggen 2010) and also physical energy component, 
such as magnetic field and relativistic particles, should be impor- 
tant to attach this problem (e.g. Bruggen & Kaiser 2001; De Young 
2003; Brighenti & Mathews 2003; Ruskowski et al.2007; Sijacki 
et al.2008; Xu et al.2008; O'Neill & Jones 2010; Scannapieco & 
Bruggen 2010; De Young 2010). 

4.2.3 Hybrid external and internal extra-heating models. 

Early (z ~ 3 - 10) pre-heating and late (z < 2) AGN feedback 
models account for a variety of energy exchanges between active 
galaxies and the diffuse baryon gas finally forming a galaxy cluster. 
The main physical difference between the two regimes is that pre- 
heating acts as an external heating mechanism, modifying the en- 
tropy of baryons in a pre-collapse phase, while AGN feedback is an 
internal heating mechanism acting within the already formed DM 
potential well of a massive halo. The two scenarios imply a very 
different energetic budget, since for the same given entropy level a 
larger energy per particle is required at higher cosmic density (e.g. 
Tozzi & Norman 2001; Mc Carthy et al.2008). Only quite recently 
cosmological numerical simulations have achived sufficient reso- 
lution and complexity to follow the interplays between the ICM 
and the populations of galaxies in a self-consistent way, along the 
whole cosmic evolution (e.g. Sijacki et al.2008; Teyssier et al.2010; 
Dubois et al.2010; Mc Carthy et al.2009). 

In order to match the two approaches in the same cluster run, 
we investigated a re- simulation adopting the intermediate (So = 
lOOkeVcm 2 , run PH2) scenario for early pre-heating, and less 
powerful jets at z = 1, with 6j e t ~ 2 • 10 57 ergs (run B5). 

In Fig|20]we show how the inner gas density and gas entropy 
of run PH2 (top panel) are modified by the late jets activity (bottom 
panel): a vertical structure of gas with entropy S ~ lOOkeVcm 2 
is found at the opposite sides of the cluster center, resulting in sig- 
nificant entrainment of the cold and dense gas of the cluster core, 
which is uplifted to larger radii at z = 0. In consequence of this, the 
inner density core is significantly depleted compared to run PH2 at 
the same redshift. 

In Fig|2T]we compare the profiles of gas density and gas tem- 
perature for the non-radiative run (R0, dashed black line), of the 
radiative run with strong pre-heating (run PH4, in red) and of run 
B5 (blue). Since the moderate amount of early pre-heating already 
prevented over-cooling of the gas in the cluster core for z > 1, the 
late injection of jets do not drive of shocks stronger than M > 2, 
and the extra entropy input is more uniformly released within the 
cluster core through mixing, rather than through violent shock heat- 
ing. Indeed, the gas entropy profile and gas density profile are very 
similar to the non-radiative case (run R0, solid blue line), and a 
well defined entropy floor is recovered; the similarity is even more 
evident at z = 0.3 (dashed blue line), ~ 2Gyr after the end of the 
jet injection. 

We conclude that, when the thermal properties of the cluster 
are concerned, the same configuration produced in a non-radiative 
run can be approximately achieved also in a radiative simulation, 
thanks to the combined effect of a uniform external pre-heating of 
So ~ lOOkeVcm 2 at high (z ~ 10) redshift and of a later phase of 
internal heating from jets injected by AGN at z ~ 1, characterized 
by an average power of Wj e t ~ 2 • 10 41 erg/s. 

This is in line with results obtained with semi-analytical 1- 
D calculation presented in Mc Carthy et al.(2008), and illustrates 
one of the (likely many) possible combinations of external and in- 
ternal heating mechanisms and radiative cooling in realistic galaxy 
cluster simulations. Of course this result is derived only from a sin- 



gle cluster object only, and additional re- simulations considering a 
wide range of masses and dynamical state are necessary to investi- 
gate the above issues even in a statistical sense. 

We note however that the kinematic structures of the ICM ve- 
locity field in the various run are very different (Fig|22l): the non- 
radiative run RO is characterized by a peaked velocity profile, with 
an infall velocity of ~ — 200km /s within the cluster core, due 
to the presence of a crossing satellite. On the other hand the other 
runs with additional heating and cooling do not present this feature; 
this is due to the fact that the same gas clump have been destroyed 
in the past by the action of early-preheating and enhanced shock 
heating while crossing the main cluster virial radius. In addition, 
the run with moderate jet feedback (run B5) shows a sharp veloc- 
ity structure in the profile of radial velocity with +250/cra/s at 
~ 200 — 400kpc/h, and a quite flat inner velocity profile, after the 
forcing of AGN outflows in the past. We note that similar features 
in the radial velocity field, in response to AGN feedback, have been 
recently reported by Dubois et al.(2010) for cosmological simula- 
tions with the RAMSES code. 

This kinematic differences should lead to different large scale 
patter of mixing/metallicity, and could be used to discriminate 
among degenerate thermodynamical structures of galaxy clusters, 
to be compared with real observations provided by future high res- 
olution spectroscopic observations of the ICM. 



5 DISCUSSION AND CONCLUSION 

In this paper we presented a detailed numerical study on the nu- 
merical and physical reasons for the presence of a flat entropy core 
in the innermost region of galaxy clusters, simulated in cosmolog- 
ical grid simulations adopting AMR (with customized version of 
the code ENZO). To this end, we performed 30 cosmological re- 
simulations of the same non-major merger cluster of final mass 
M « 3.1 • 10 14 M©//i. We accurately studied the many parame- 
ters likely affecting the production of physical or numerical entropy 
within cosmological cluster simulations. 

Evidences were presented for the ubiquitous presence of a 
well defined entropy floor inside the cluster core radius of non- 
radiative cosmological runs, mostly independently on the numeri- 
cal details of the simulations (Sec(3}. This plateau-like entropy dis- 
tribution ranges from ~ 70keV • cm 2 to ~ 150keV • era 2 , and has 
a size of r CO re ~ 100 — 200kpc/h (~ 0.1R V i r ) at z — 0. The most 
relevant factor affecting the level of entropy in the cluster core, 
among those investigated, is the mesh refinement strategy adopted 
in the simulation (Se d3.21 ): when mesh refinement triggered by ve- 
locity jumps is added to the standard mesh refinement triggered by 
gas/DM over-density, the level of entropy in the core is increased 
by a factor ~ 1.5 — 2. This is due to the enhanced presence of 
mixing motions and shocks within the cluster, which are otherwise 
more dampened by under- sampling effects in a standard refinement 
strategy. The other effects (N-body gravitational noise, cold unre- 
solved flows, softening length in the calculation of DM dynamics, 
spatial resolution for gas dynamics) were found to introduce dif- 
ferences at the order of a few tens of percent on the final entropy 
profile, without changing the inner slope of the radial distribution 
in a relevant way (Sec l3.31 Sec l3.51 ). 

In Sec l4.ll we explored in detail the physical mechanism 
which produces the flat entropy distribution in non radiative cosmo- 
logical simulations, using Lagrangian tracers advected in the sim- 
ulation. The regular accretion of shells of matter onto the forming 
cluster is responsible for the strong entropy stratification found for 
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Figure 20. Gas density (colors) and gas entropy (contours spaced in 
Alog(S) = 0.2) for central region of run PH2 (top panel) and run B5 
(bottom panel). The side of the images is 2.5Mpc/h; the colors are as in 
FigEl 

r > lOOkpc/h at the final epoch; this stratification sets up already 
in the first phases of the collapse (z ~ 1 — 2), and mirrors the differ- 
ent thermodynamic history of clumpy and smooth accretions onto 
the main cluster. Inside r < lOOkpc/h, mixing motions driven by 
matter accretion gradually mix lower and higher entropy gas, pro- 
ducing an almost constant entropy profile in the region where the 
gas+DM gravitational potential is shallow. Our results here con- 
firm and extend the seminal work of Mitchell et al.(2008) to a fully 
cosmological framework, and to non-major merger galaxy clusters. 

We also explored more complex physical modeling of the 
ICM, following the effect of radiative cooling (Sec l4.2l ) and of 
non-gravitational heating mechanisms, such as early uniform pre- 
heating (Sec l4.2.l1 ) or late jets-like injection by AGN (Sec l4.2.2l ). 
We report that, while it is feasible to recover a very similar profile 
of non-radiative runs and of observed CHANDRA clusters (Cav- 
agnolo et al.2009) with a suitable choice of uniform extra-entropy 
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Figure 21. Gas density and entropy radial profiles for the run PH2 (red), 
run B5 at z = (solid blue) and z = 0.3 (dashed blue). The additional 
black lines shows the profile for the non-radiative run (run R0, long dashed) 
and for the cooling run (run R15, solid). The vertical dotted arrow shows 
the softening length adopted in run R0, while the dot-dashed one shows the 
softening of runs PH4 and B5. 

input at z ~ 10 (in the range of So ~ 100— 200keVcm 2 ), it results 
impossible to achieve similar results with the set of trials of only 
jet-like injections, where we simulated the thermal feedback from 
a central AGN. In the explored configuration, the main problem 
is that the bulk of the energy release from AGN outflows triggers 
too intense shock heating in the cooling ICM at z ~ 1, and steep- 
ens the inner gas entropy profile as shown in semi- analytical model 
presented in Lapi, Cavaliere & Menci (2005). However, the adop- 
tion of a hybrid model with moderate early pre-heating and late and 
modest jets feedback (Sec l4.2.3l) is found suitable to recover a ther- 
modynamic structure which extremely similar to the non-radiative 
fiducial run (R0), and within observations (Cavagonolo et al.2009). 
This suggests that the cluster configurations generally produced in 
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Figure 22. Profiles of the total radial (solid lines) and chaotic radial (dashed 
lines) module of velocity for run R0 (red), for run PH4 (blue) and for run 
B5 (green). The dotted arrow shows the softening length adopted in run R0, 
while the dot-dashed one shows the softening of runs PH4 and B5. 



cosmological non-radiative runs may be considered, at first approx- 
imation, a viable representation of real galaxy clusters with cooling 
and feedback mechanisms at work. However, similar thermal dis- 
tributions at z = may be characterized by quite different kine- 
matic structure, depending on the different feedback mechanisms 
at work, leading to potentially detectable spectroscopic X-ray fea- 
tures. The study of other important observables not considered in 
this paper (such as the distribution of stars, metals, and the cluster 
baryon fraction) is expected to provide additional ways to discrim- 
inate between similar thermal models (e.g. Booth & Schaye 2009; 
Teyssier et al.2010). 

Our conclusions is presently limited by the fact that we 
adopted so far "ad-hoc" models of pre-heating and AGN feedback, 
which are coupled to the simulation "by hand", and that we focused 
on only two clusters (see also the Appendix). It would be interest- 
ing to investigate the same issues using a number of clusters with 
different masses/dynamical history. 

As a final remark, we note that the mechanism which produces 
the entropy floor in cosmological grid simulations of non-radiative 
clusters can also explain the long-debated difference reported for 
SPH and grid run (e.g. Frenk et al.1999; O'Shea et al.2005; Tasker 
et al.2008; Wadsley et al.2008; Mitchell et al.2008; Springel 2010). 
In major merger clusters, the difference in the core entropy is set by 
different efficiency in the mixing of gas at the moment of the closest 
encounter between the clusters, as convincingly shown by Mitchell 
et al.(2008). In relaxed clusters, a similar mechanism works on 
longer time scales, due to the continuous action of subsonic chaotic 
motions triggered by the accretion of satellites. Any different in 
the modeling of mixing in the two numerical methods can explain 
the presence or absence of a well defined entropy core structure. It 
has been shown that the presence of an artificial viscosity term in 
standard SPH greatly reduces the small scale mixing in a number 
of realistic cases, compared to standard Eulerian simulations (e.g. 
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Agertz et al.2007; Wadsley et al.2008; Springel 2010), and that the 
adoption of less viscous simulations produces entropy distributions 
in clusters more similar to Eulerian runs (e.g. Dolag et al.2005; 
Mitchell et al.2009). As shown in this work, shock heating is the 
leading source of entropy production in cluster (well beyond the 
role of any possible numerical artifact), while physical mixing is 
the reason for the spreading of entropy in the innermost cluster re- 
gions. The fact that mixing in SPH is usually reduced by numerical 
effects, fully explains while the two methods are in disagreement 
in the center of clusters, when mixing is maximum in grid codes, 
while they are found in much better agreement at larger radii (e.g. 
Frenk etal.1999). 

Since the buoyancy in the stratified ICM is strongly depen- 
dent on the underlying entropy distribution, the above findings em- 
phasize the need of having a suitable numerical representation of 
cluster cores, since this may affect also the estimated energy bud- 
get needed from non-gravitational heating mechanism, and in their 
efficiency in mixing/heating the surrounding Intra Cluster Medium. 
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7 APPENDIX 

Cluster mergers may boost shock heating and mixing motions in 
the ICM for several Gyrs (e.g. Ricker & Sarazin 2001), significanly 
changing the physical entropy generation in a cluster with a dynam- 
ical evolution different than the relaxed case explored in the main 
part of the paper. 

We present here some complementary tests on the entropy 
distribution of a major merger cluster of final mass M « 2.1 • 
10 14 M©//i. The bulk of the total mass of this cluster is assembled 
in a major merger at z = 0.85, with an approximate mass ratio 
of M1/M2 ~ 3 between the colliding halos. Only a subsample 
of the re- simulations presented in the main part of the paper were 
repeated with this cluster; the parameters of the tests run in this 
case are listed in Table [7] (all the cosmological parameter are as in 
Sec l3.ll ). The non radiative simulation of this cluster produces a flat 
entropy core with S ~ 80 — 90keV cm 2 forr < lOOkpc/h, similar 
to the non-radiative fiducial run (run R0, Sec l3.21 ). 

The adoption of the additional mesh refinement based on ve- 
locity jumps (run A2) causes a net increase of the internal entropy 
compared with the more standard mesh refinement strategy based 
on gas/DM over-density (Al), as shown in Sec l3.2l The role of the 
gravitational softening is found to be more important in this major 
merger cluster, compared to the relaxed cluster studied in the main 
part of the paper: the adoption of a larger softening (50kpc/h, run 
A3) produces an entropy core larger by ~ 20 per cent compared to 
a smaller softening (25kpc/h, run A2). This stresses the higher im- 
portance of having a good resolution for the computation of gravi- 
tational forces in the case of violent oscillations of the gravitational 
potential driven in a merger event, which may generate an amount 
of extra-entropy production of numerical origin. 

The adoption of radiative cooling (in the case of standard re- 
finement, CI, or with the velocity-based refinement, C2) causes a 
very similar trend as in the case of the relaxed cluster explored in 
the paper, with the onset of catastrophic cooling for r < 200kpc/h. 
This shows that, at least for this early major merger (z ~ 0.85), 
the action of intense heating from merger shocks is not effective in 
destroying the forming cooling region, somewhat at variance with 
other works with Eulerian simulations (e.g. Burns et al.2008). 

To spare computational time, most of the runs with non- 
gravitational heating were performed only with the standard re- 
finement scheme. A uniform pre-heating of So — lOOkeVcm 2 at 
z = 10 (run PH2) is ineffective to stop the catastrophic cooling for 
r < 50kpc/h; however when applied to this major merger system 
it results in a significantly higher inner entropy value compared to 
the relaxed system studied in the main part of the paper (Sec l4.2.l1 ). 
A pre-heating of So = 200keVcm 2 (PH3) on the other hand al- 
most perfectly recover the entropy distribution of the non-radiative 
run (PH3) for the merger cluster. 

Runs with pre- heating of S = lOOkeVcm 2 and AGN feed- 
back were produced for the jets energy of 6j e t = 2 • 10 57 ergs 
(run Jl) and 6j e t = 10 58 ergs. The profile of the Jl run is very 
similar to what obtained for the relaxed cluster of the paper; the 
profile with a higher AGN energy results in a flat entropy profile at 
~ 130keVcm 2 for r < 200kpc/h. Finally, we re-simulated run 
J2 adopting the velocity based refinement (run J3), finding still a 
very flat profile inside r < 200kpc/h, and a ~ 50 per cent larger 
entropy in the center. According to the assumed jet energy, our re- 
simulations with cooling, pre-heating and AGN feedback can thus 
provide an acceptable match with one of the two classes of the bi- 
modal entropy distribution reported for the CHANDRA observa- 
tions of Cavagnolo et al.(2009). 
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Table 2. Main characteristics of the performed runs (as in Tab|2j.In the last 
row, the assumed pre-heating background is So = lOOkeVcm 2 and the 
thermal energy of the jets is 2 • 10 57 ergs 

ID Max Res. [kpc/h] soft, [kpc/h] AMR note 



Al 


25 


25 
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non-radiative 


A2 


25 


25 


DV 


non-radiative 


A3 
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50 


DV 


non-radiative 


CI 


25 


25 
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cooling 


C2 


25 


25 


DV 


cooling 


PH2 


25 


25 
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cool.+PH(100fceVcm 2 ) 


PH3 


25 


25 
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cool.+PH(200/eeV cm 2 ) 


Jl 


25 


25 
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cool.+PH+J(2 • 10 57 ergs) 


J2 


25 


25 
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cool.+PH+J(10 58 ergs) 


J3 


25 


25 


DV 


cool.+PH+J(10 58 ergs) 




Figure 23. Profiles of the gas entropy in the major merger run. The ad- 
ditional long dashed line show the bimodal gas entropy distribution from 
CHANDRA observations (Cavagnolo et al.2009). 



In Fig[7]we show the time evolution for three snapshots of run 
CI, J2 and J3, showing the X-ray bolometric luminosity for a region 
of ~ 1.8 x 2.2Mpc/h per side and thickness 25kpc/h, around the 
epoch of the major merger. Comparing CI to J2, we show how 
the action of AGN feedback removes the central gas condensation 
within the cooling region in a few ~ lOMyr. Run J3 emphasizes 
the role played by the mesh refinement strategy on the expanding 
shocks driven by the (almost contemporary) AGN burst and the 
major merger. We also note how at least one inflated "bubble" can 
survive for a few time steps after the injection, due to the reduced 
numerical mixing in the implemented mesh refinement strategy. 

The above tests suggest that the most important findings re- 
ported in the main body of the article are general, since they apply 
to clusters with a similar mass but two completely different dy- 
namical histories. However the efficiency of the extra-heating mod- 
els applied to radiative runs may depend on the dynamical history 
of the host cluster, and further studies are needed to estimate the 
global efficiency of the proposed scenarios in a statistical sense. 
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Figure 24. Maps for the evolution of the X-ray bolometric luminosity for slices of 1.8 x 2.2Mpc/h and depth 25kpc/h for the major merger cluster studied 
in the Appendix, at z = 0.9, z = 0.85 and z = 0.81. The top panels are for the pure cooling run (CI), the middle panels are run Jl with cooling, pre-heating 
(S = lOOkeVcm 2 at z = 10) and AGN feedback (e je t = 10 58 erg/s), while the bottom panel are for a re-simulation with the same setup, but mesh 
refinement triggered by velocity jump (run J3). 



